Land use and land cover changes and their impact on ecosystem service values in the north-eastern highlands of Ethiopia

The land use and land cover (LULC) changes driven by the growing demands of mankind have a considerable effect on ecosystem services and functions. The study was carried out in the north-eastern highlands of Ethiopia to (1) analyze the effect of LULC changes between 1984 and 2021 and (2) assess the spatiotemporal variations in ecosystem service values (ESVs) and elasticity in response to LULC changes. Using Landsat imageries from 1984 to 2021, the spatiotemporal changes in LULC were evaluated with supervised image classification using maximum likelihood algorithm in ArcGIS software. Six LULC types were subsequently categorized, with overall accuracy and Kappa coefficients above 87% and 0.87, respectively. The ESVs were then estimated based on the Benefit Value Transfer (BVT) approach employing modified conservative value coefficients. The findings revealed a significant increase in cultivated land (9759.1ha) and built-up area (10174.41ha) during the stipulated periods and a drop in other land use types. The forest loss gradually decreased from 4.1% in the second period (1991–2001) to 0.58% in the third (2001–2021), compared to the first of the 1.1% conversion rates. Similarly, the proportion of grassland and water bodies steadily reduced over the stipulated periods, by 1.15% and 2.3% per annum, respectively. The overall loss of ESVs in the study landscape was estimated to be 54.4 million US$ (67.3%), drastically decreasing from 80.3 million US$ in 1984 to 26.4 million US$ in 2021, driven by the declining area coverage of water bodies, grassland, and forestland. Regardless of the loss, the ecosystem functions of hydrological regulation (37.2, 35.0, 6.1, and 5.1 US$ ha-1yr-1), water supply (14.5, 13.6, 2.4, and 2 US$ ha-1yr-1), and food production (9.8, 10.0, 9.1, and 9.9 US$ ha-1yr-1) contributed the most to the total ESV of each year while disturbance regulation and cultural values were the least throughout the study periods. The coefficient of sensitivity (CS) analysis revealed that our estimates were relatively robust. The findings further showed that human-dominated land-uses at the expense of natural ecosystems are the primary drivers of LULC transitions and the ensuing loss of ecosystem services in the region. Thus, this calls for intensive work on more effective land use policies that encourage an integrated management approach, with a focus on safeguarding the sustainability of natural ecosystems.


Introduction
Tropical forest ecosystems provide a wide range of services including habitat conservation, water and soil conservation, biodiversity preservation, and carbon sequestration [1,2].However, nearly half of the tropical forest that existed at the beginning of the century has already vanished [3].Land-use change is thought to be the major driver of the loss of forest ecosystems and the change in the accompanying services [1,4].Changes in land use/land cover (LULC) have a marked effect on local-to-global ecological processes and ecosystem services [2,4], pose a serious threat to biological resources, such as forests, grasslands, wetlands, and other ecosystems [1,5], and have a profound impact on achieving Sustainable Development Goals (SDGs) [5,6].Collectively, all of which have a strong influence on carbon sources and sinks as well as global environmental change.The main driving forces of the LULC transition might be connected to multiple factors, such as the increasing human population and accompanying consumption demands [7,8], rapid socio-economic development [9,10], policies and institutional factors [10,11].Climate-induced changes (e.g., drought, precipitation variability, forest fire, invasive species, and other disturbances) can also put pressure on LULC alterations [2,12,13].
In Ethiopia, LULC changes especially forest cover transition has been more severe in the highland areas (altitude >1500 m) over the past few decades [14][15][16][17][18][19][20][21], as a result of increasing population, which led to significant forest clearance for agricultural use, overgrazing, etc.The forest under this threat covers more than half of the Afromontane forests of the country [14], which is one of the key biodiversity hotspots in Eastern Africa with a diverse flora and significant ecological reservoirs [12,22].According to FAO [23], Ethiopia's forest cover decreased from 15.11 million ha in 1990 to 12.5 million ha in 2015.Several recent studies have also shown that the proportion of forested land reduced by 8.2% in the Northern [24], 3.5% in the North-western [17], 34.45% in the Western [25], over 45% in the south western [26], and 13% in the central highlands of Ethiopia [21].A substantial decline in forest cover has been attributed to its conversion to abandoned pasture and croplands.In most cases, immense pressures from anthropogenic activities (e.g., poor agricultural practices, deforestation, overgrazing, etc.) coupled with a rapidly increasing population are the main drivers for changes [14-18, 24, 25, 27].As a result, forest resources, mainly natural forests, are increasingly becoming scarce [14,22], resulting in habitat fragmentation, biodiversity loss, and a reduction in ecosystem services.
Ecosystem services refer to the benefits that humans derive from natural ecosystems and their components, including provisioning services (e.g., food, timber, fodder, water, medicinal products), regulating services (e.g., carbon storage and sequestration, water purification, preventing soil erosion), supporting services (e.g., habitat, biodiversity conservation), and cultural services (e.g., aesthetic, recreational, and tourism) [5][6][7][8].Human well-being and the functioning of the global economy depend on ecosystem services, but these services are under threat because of the intricate interplays between people and the environment, which result in ecosystem degradation and biodiversity loss [28], and lead to economic impoverishment.It is estimated that 60% of ecosystem services have changed in just the last 50 years [5,6], with humancaused LULC alterations being the primary causes [5][6][7][8].Even though ESVs vary spatially, Costanza et al. [8] estimated that global LULC changes have caused ESVs to decrease from US $ 145 trillion year −1 in 2007 to US$ 125 trillion year −1 in 2011.The decline in ESVs is more pronounced and burgeoning in developing countries like Ethiopia [29][30][31].Recent studies in Ethiopia [29,31,32] have asserted that the changes in ESVs resulting from LULC changes have led to a significant reduction in forest coverage and an increase in croplands.According to Tolessa et al. [30], forest conversion has resulted in a loss of 40.7% (3.69 million) in ESVs in the central highlands of Ethiopia, and comparable results have been found in other parts of the country [32,33].
To avert the scenarios of forest degradation and ESVs loss, Ethiopia's government has adopted the Green Legacy Initiative through sustainable forest management and restoration programs, which aims to address the issue of forest degradation and loss of ecosystem services through sustainable forest management and restoration programs [34].These initiatives are designed to monitor and manage regional ecosystem changes, provide a framework for ecological restoration and conservation, and also form sustainable development policies, with a focus on the highlands [34].The North-eastern highlands of Ethiopia are notable for having mosaics of different LULC types where dry evergreen Afromontane forests have historically been degraded as a result of human pressure.However, there have been recent efforts by the government to restore these forests, including community-driven initiatives such as tree planting, enclosures, and exclosures, as well as soil and water conservation measures (SWC) for over two decades.Despite these interventions, previous research has not investigated the impacts of LULC changes on ESVs in the region, considering the differences in landscapes and restoration endeavours.Therefore, this study was initiated to analyse the (1) trends in LULC changes/trajectories and vegetation covers of the area from 1984 to 2021, (2) spatiotemporal variation in ESVs of the area caused by LULC change in the stipulated periods, and (3) contribution of individual ecosystem service values and elasticity driven by landscape dynamics during the study periods.The findings of the study were instrumental to show the environmental and ecological changes, project additional intervention options and provide evidence-based decisions concerning regional sustainable development.

Description of the study area
The study area falls into two districts, namely Guba Lafto and Woldia administrative city.The areas are part of the North Wollo administrative zone of the Amhara National Regional State (ANRS), in the Northern eastern Ethiopian highlands (Fig 1 ), and are located between 112 0'00''N to 12˚10'.00"Nlatitude and 39˚20'00'' to 39˚40 ' 40.00"E longitude.The study area is part of the north-central plateau, consisting of a chain of mountains, hills, valleys, plains, and plateaus with topography that ranges from nearly 1379 to 3200 meters above sea level (m.a.s.l) [35].The study area has predominantly midland (Woina Dega) and highland (Dega) agro-ecological zones, with a bimodal rainfall pattern (from February to March and from June to August).It has an average annual rainfall ranging from 10 mm in January to 380 mm in August and temperatures of 12.4 to 28.7˚C [36,37].The major soil types of the area are the vertisols, nitosols, and cambisols [38].The main economic activities are farming and livestock rearing, and the dominant crops are Sorghum bicolor, Zea mays, Eragrostis tef, and Cicer arietinum, but the area is often highly degraded.Despite centuries of human use and land modification, the natural vegetation of the study area belongs to the remnant dry evergreen Afromontane forests of the Ethiopian highlands [39].Natural forests can be found on steep slopes and in churchyards, dominated by Olea europaea subsp., cuspidata, Eucalyptus camaldulensis, Acacia spp, etc.

Data acquisition
Changes in land use/land cover (LULC) are the major factors affecting ecosystem processes and services [29].The changes are happening at a fast rate due to an increase in human demands on the forest ecosystem [40].Thus, the evaluation of long-term regional and global LULC change is made by the Landsat data archives, which are imperative for historical records of the Earth's land surface [41].For this study, multispectral satellite imageries were downloaded from the United States Geological Survey Earth Explorer website (USGS) (https:// earthexplorer.usgs.gov/) to analyse LULC changes in the study area over 37 years (Tables 1  and 2).Hence, Landsat 4-5 Thematic Mapper (TM) of 1984, Landsat5 Thematic Mapper (TM) of 1991, Landsat7 Enhanced Thematic Mapper Plus (ETM+) of 2001, and the Landsat8 Operational Land Imager (OLI) of 2021 were used for this study (Table 1).The study was purposely commenced in 1984, as it aligned with a time when deforestation was rampant, particularly in the dry evergreen Afromontane regions of the country [42].It is crucial to ascertain the  magnitude of the issue in the region and devise a proactive decision-making process towards solving intricate environmental problems.The images were captured with high-resolution (30m x 30m) during cloud-free scenes of the dry season (December to January) to avoid the effects of seasonal variations in ground cover especially effects of annual crops or cloud cover, which might affect the data interpretation [17,41,43], and are particularly valuable for studying natural and agricultural environments, as they can enable informed decision-making that promotes sustainable land use and healthy plant communities [17,43].These images were projected to the Universal Transverse Mercator (UTM) using the World Geodetic System (WGS) 84 zone 37 North datum.In addition, auxiliary data (ground truth points) were acquired between September and October of 2022 using a portable handheld Global Positioning System (GPS) and a Google Earth satellite image for the study area.Hence, a total of 350 GPS points (50 points per class) were used to classify and validate LULC classifications.

Image pre-processing and analytical methods
All imageries underwent pre-processing (e.g., layer stacking, mosaic, geometric and radiometric calibration, and panchromatic) using ArcGIS (version 10.7.1 software, Esri, California, USA) to enhance the image quality.The radiometric calibration was carried out using the raster calculator (Spatial Analyst) tool.However, the geometric adjustment was not required because the pixels in the Landsat images from the USGS website are accurate [44].

Normalized Difference Vegetation Index (NDVI)
The Normalized Difference Vegetation Index (NDVI) is one of the most often used metrics for detecting vegetation conditions (e.g., cover, biomass, primary production, and carbon balance) [44][45][46].It also allows for the evaluation of spatial and temporal changes occurring on the Earth's surface using reflectance measurements [45,47].In this study, NDVI values were calculated based on the reflectance of the red (R) and near-infrared (NIR) bands (Eq 1) of satellite images (https://earthexplorer.usgs.gov/).

Supervised images classification
Supervised image classification was employed using the maximum likelihood algorithm, which is a widely used technique in image classification [49].Ground truth points from each LULC category were collected to aid in the classification process.The supervised classification method is commonly used to identify land use classifications using spectral signatures [50], and it assumes that the statistics are normally distributed for each class and in each band [51].
The image was then sorted into groups based on how closely it resembled the training signatures after collecting samples from the various land cover classes (training data) [52].The variance and covariance of the spectral response patterns were quantitatively evaluated using the maximum likelihood algorithm, and each pixel was then assigned to the class for which it has the greatest likelihood of association [41].The LULC classes that were identified are forestland, cultivated land, grassland, built-up areas, barren land, and water bodies (Table 2).

Detection of land use/land cover changes
A post-classification approach was utilized for identifying LULC changes over time based on the numerical values extracted from classified images [41,44,51,53].This technique involves overlay comparisons to determine which LULC types have been converted [41,54].To do this, classification images of the same region generated independently at distinct time periods are compared.Then, the final classified thematic maps from subsequent periods (i.e., 1984-1991, 1991-2001, 2001-2021, and 1984-2021) were cross-tabulated [29,33,41,44] to determine the percentage and rate of land use/land cover change detections over time series.To further understand the observed LULC transitions, the total and net changes of the area were calculated using Eq 2.
Where, LULC I is the degree to which LULC changed or transformed during the reference period i, LULC end year and LULC start year represent the area of a particular land use unit at their respective beginning and ending times, and t denotes the estimation period.Thus, positive values indicate an increase in the extent of LULC, whereas negative values show a decrease.

Accuracy assessment
The accuracy assessment indicates the degree to which a particular land cover on the ground deviates from the reference map or the frequency with which actual features on the ground are accurately depicted on the classified map [41,55].All of the classified imageries in the current study had their accuracy evaluated using ground truth points collected from varied land use types in the eco-region via: (1) direct data acquired using GPS (Garmin GPSMAP 60CSx) from the field, (2) information obtained from interviews with elders and communities in the study area (3) Woody Biomass Inventory and Strategic Planning Project data (https://www.acronymfinder.com/Woody-Biomass-Inventory-and-Strategic-Planning-Project-(Ethiopia)-(WBISPP).html), and (4) Google Earth data sets.As a result, the reliability of the classified LULC was assessed using confusion matrices: the user's accuracy, producer's accuracy, overall accuracy (Eq 3) [51] and kappa coefficient (Eq 4) [56] as follows

Overall accuracy
Where, Xii is the diagonal elements in the confusion matrix, X is the total number of samples in confusion matrix Where, K = kappa statistics, M = total number of observations in the matrix, r = number of rows in the confusion matrix, n ij = number of observations in row i, column j, n i = total number of observations in row I, n j = total number of observations in column j.According to Viera and Garrett [56], values of kappa statistics range from -1 to 1, i.e., values greater than 0.80 denote substantial or strong agreement, values between 0.40 and 0.80 represent moderate agreement, and less than 0.40 signify poor agreement between classification and reference data.

Estimation of ecosystem service values (ESVs)
Ecosystem service (ES) valuation signifies the monetary value assigned to a variety of ecosystem services offered by the natural environment [8,57,58].This value is predicated on the assumption that the average unit value across all locations is homogeneous, which may not be the case in reality [8,29,58].In this study, the benefit value transfer (BVT) approach was employed to estimate the ecosystem service values (ESVs) in response to the spatiotemporal LULC change of different land uses for the years 1984, 1991, 2001, and 2021.The BVT allows estimates for one site to be based on estimates from other sites or data from studies that have previously been conducted in different places [29-33, 59, 60].Costanza et al. [56] have classified the whole biosphere into 16 biomes and 17 service function types (Tables 3 and 4).The ESVs of each LULC types at various locations are estimated, and thus, the equivalent value per unit area of ES for the study area was extracted using the modified conservative value coefficients, based on data from Kindu et al. [29], which was modified after Costanza et al. [8,57] models.This is because the study area has comparable similar market values and geographical settings, notably tropical Afromontane region (e.g., 1500-3000 m.a.s.l).Accordingly, ESV per unit area for each LULC category (Table 2) was assigned based on the corresponding coefficient value of the ES (Tables 3 and 4).Land use classifications such as barren land and settlement did not have a coefficient in other studies [e.g., 8, 29-33, 57, 58].Similarly, no coefficients were assigned to built-up areas and barren land in this study since the LULC classes lacked vegetation.The ecosystem service value (Eq 5), ecosystem function (Eq 6), and total ESV (Eq 7) for each thematic class were computed after evaluating the ESV per unit area for each land cover class using the models described by Costanza et al. [8,57] and recommended by Kreuter et al. [61] and Kindu et al. [29].
Where, ESV k , ESV f , and ESV refer to the ecosystem service value of LULC type 'k', service function f, and the total ESVs, respectively (in unit of US$).A k represents the area (ha) for LULC category 'k' and VC kf is the value coefficient for land use type k with ecological service function type f.
The changes in ESVs were evaluated by calculating the difference between the estimated values in each reference year (Eq 8).

Analysis of coefficient of sensitivity
The sensitivity analysis was used to verify the ESV variations in response to ±50% adjustments of the ESV coefficients for each type of LULC [29,62].The coefficient of sensitivity (CS) was therefore determined based on the standard economic concept of elasticity of formula noted by Kreuter et al. [61].It follows as: Where, ESV i and ESV j are the initial and adjusted total estimated ecosystem service values, respectively.VC ik and VC jk are the initial and adjusted value coefficients (US$ ha −1 year −1 ) for LULC type 'k'.The estimated ecosystem service value is, however, elastic for the value coefficient if the changes of the two values are greater than the threshold limit (> 1), and it is fully inelastic when the value is less than one [61].

Accuracy assessment of LULC maps
Table 5 displays the accuracy assessment for the identified LULC types.For the years 1984, 1991, 2001, and 2021, the overall categorization accuracy ranges from 87 to 91%.According to Monserud and Leemans [63], a very good kappa value should fall between 0.70 and 0.85.In this study, the classified images exhibited a kappa coefficient above 0.80.As a result, the validation data set revealed a very strong agreement between the categorized images and the ground truth data.

Land use land cover (LULC) changes
Trends of LULC trajectories, 1984 to 2021.The LULC analysis revealed that the study areas had noteworthy changes in land use and cover during the study periods .The LULC change trajectories in Figs 2 and 3 represent an increasing trend in the proportion of cultivated land, built-up area, and barren land, which make up roughly 78.9%.However, the other three LULC types (forestland, grassland, and water bodies) had decreased during the study periods, together comprising 21.1% of the total (Figs 2 and 3).Notably, the proportion of forestland had decreased over the past 37 years; from 10.8% (10706.4ha) in 1984 to 9.7% (9631.84ha), 5.63% (5583.04ha), and 5.0% (5005.31ha) in 1991, 2001 and 2021, respectively, with 1.44% (-154.1 ha) of an annual rate of deforestation (Figs 2 and 3).Similar trends were also seen in the grassland, which decreased from 26.5% (26319.6 ha) in 1984 to 24.6% (24451.17ha), 19.13% (18985.62ha), and 15.2% (15051.85ha) in 1991, 2001, and 2021, respectively, with an annual loss of 1.2% (Figs 2 and 3).Additionally, the area covered by water bodies declined drastically from 6.9% (6806.7 ha) in 1984 to 914.89 ha (0.9%) in 2021, at a rate of 2.34% (159.24ha) per annum (Fig 3).The analysis of LULC changes showed the most significant land cover instability over the last 37 years.During the first period (1984-1991), there were changes observed in the LULC.The areas used for cultivated land, built-up, and barren land increased by 0.89% (2093.53ha), 0.63% (1074.56ha), and 0.87% (833.93 ha), respectively, indicating annual rates of change of 0.13%, 0.09%, and 0.12% (Figs 3 and 4).Conversely, the percentage of forest cover, grasslands, and water bodies declined during the same time frame, with an annual rate of change of -0.2% (-153.5 ha) for forestland, -0.14% (-266.9ha) for grassland, and -0.12% (-57.3 ha) for water bodies, as depicted in Figs 3 and 4. From 1991 to 2001(second period), there was an increase in the area covered by cultivated land (384.48 ha; 0.11%), built-up area (7985.18ha; 8.24%), and barren land (6439.77ha; 4.5%), at a rate of 0.011%, 0.82%, and 0.45% per annum, respectively (Figs 3 and 4).At the same time frame, the percentage of forest cover, grasslands, and water bodies showed a declining trend, where the annual rates of change was -0.42% (-404.9ha), -0.22% (-545.6 ha), and -0.83% (-529.5 ha), respectively (as used in Figs 3 and 4).Between 2001 and 2021 (third period), the proportion of cultivated and built-up area showed an incremental trend compared to the transformation of the other three LULC types that exhibited the same declining pattern as in the first and second periods (Figs 3 and 4).In contrast to the changes observed between 1991 and 2001, barren land area decreased at a rate of -1.05% (-4346.56ha) (Fig 3).Additionally, the extent of forest area conversion was lower in the third period (2001-2021) than the first and second periods, while grassland and water bodies continued to reduce by -1.04% and -0.9%, respectively (Figs 3 and 4).The overall LULC change results (1984-2021) demonstrated the tremendous expansion of agricultural lands and settlements in the studied area, unlike the declining trend recorded for forestland, grassland, and water bodies (Fig 3).
Trends in transitions between LULC change analysis of the districts revealed that 52.98% (52,590.41ha) of the land underwent continuous land cover conversion while 47.02%  6).The LULC change matrix analysis in Table 6 indicated that the forestland was converted to grassland and cultivated land at proportions approximately equal to 31.5% and 12.2% over the study periods, with 0.85% and 0.33% annual rates of deforestation, respectively.In addition, the grassland was transformed into cultivated land (10432.6 ha) and built-up area (6997.4ha) at a declining rate of 0.3% and 0.2% per annum, respectively (Table 6).In 1984, out of the total area of water bodies (6806.7 ha), 6114.4ha(89.8%) was converted to other LULC types over the study period.Of this, 15.7% and 3.3% of the water body were converted to cultivated land and built-up area, respectively (Table 6).Overall results showed that grassland experienced the greatest area coverage loss.However, we found considerably higher conversions from other LULC types to cultivated land (22479.31ha-gains) and built-up area (18227.5 ha-gains) (Table 6).

Spatial-temporal variation in ESV
Effect of LULC change on the total estimated ESVs.The ecosystem service values were evaluated for the years 1984, 1991, 2001, and 2021.In 1984, the total ESV of the area was estimated to be $ 80.8 million using a modified conservative coefficient (Fig 7 and Table 7).However, the value coefficient varied by LULC type whereby, the water bodies had the highest contribution to the total ESV (55.1 million USD) during this time, followed by forestland (10.6 million USD), and grassland (7.7 million USD) (Table 7).The total ESV in 1991 was reported to be $76.4 million, in which the water bodies (51.9 million USD) were found to have the highest contribution, followed by forestland (9.5 million USD) and cultivated land.Our findings showed that the total ESV was significantly lower in the ensuing decade (2001) than in the prior research year by more than 35% (Table 7).As for each LULC type, water bodies (9.0 million USD) and cultivated land (7.9 million USD) generated the highest ESV in the area.Forestland and grasslands contributed the least to the total ESV, each with a 20% contribution (Table 7).In 2021, the total ESV was $26.4 million, with the highest contribution derived from cultivated land ($9.6 million) (Table 7).On the other hand, we considered the contribution of built-up area and bare land to the total ESV to be zero, as it was also in the modified conservative estimates (Tables 3 and 4).
Over the course of the first seven years (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991), the area's total ESV decreased by 5.4% ($4.4 Specifically, the ESVs of water bodies, forestland, and grassland decreased by $3.2 million, $1.1 million, and $0.5million, respectively, while the ESVs of cultivated land increased by $0.5 million (Table 7).However, during the 1991-2001 period, a markedly higher conversion of forestland and water bodies to cultivated and built-up areas was observed, resulting in a total ESV loss of 48.4 USD million with a proportional loss of 63.4%.Over the same period, reductions in ESVs for the water bodies (-42.9USD million), forestland (-4.0 USD million), and grassland (-1.6 USD million) were recorded (Table 7).In contrast to the noticeably higher ESV loss in the second period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001), the overall ESV loss between 2001 and 2021 was relatively lower (-1.6 USD million), with a proportional loss of -5.7% (Table 7).In this time frame, the total loss of ESV was mainly brought by grasslands (-21.4%) and water bodies (-17.8%).In the contrary, the total ESV of cultivated land increased by 1.7 USD million, with proportion changes of 21.5% (Table 7).Overall, from 1984-2021, the total ESV of the study landscape showed a downward trend, decreasing from 80.8 USD million to 26.4 USD million, resulting in a total decrease of 54.4 USD million (-67.3%) (Table 7).During this period, a substantial ESV loss of water bodies' from US$ 55.1 to US$ 7.4 million (-86.6%) and terrestrial ecosystem, i.e., forestland (-53.8%) and grassland (-42.9%) were observed over the study periods, while the cultivated land exhibited an increasing pattern from US$ 7.4 to US$ 9.6 million (29.7%) (Table 7).8. Our results indicated that there was a notable variation in the rate of change in the value of ESVf during the study periods.In 1984, hydrological regulation, water supply, food production, and waste treatment contributed most to the total ESV in the north-eastern highlands of Ethiopia.They accounted for 46.02% (37.2 US$ ha -1 yr -1 ), 17.9% (14.5 US$ ha -1 yr -1 ), 12.1% (9.8 US$ ha -1 yr -1 ), and 8.3%  (6.7 US$ ha -1 yr -1 ) values, respectively.However, disturbance regulation (0.07%) and cultural (0.02%) ecosystem service functions had the lowest contribution in the area (Table 8).Also, the same higher and lower ESVf contributions to the total ESV were recorded in 1991, 2001 and 2021 (Table 8).It is noteworthy that the rates of change of ESVf quite decreased from 2001 to 2021 compared to the more negative rates of change between the first (1984-1991) and second (1991 and 2001) phases (Table 8).Indeed, as a result of LULC changes, the entire ecosystem service functions showed a decreasing trend through the study periods, with the exception of food production and biological control, which lingered in change at almost the same or negligible rates across the study periods (Table 8).During the study period, the ESV function of hydrological regulation (-86.3%),water supply (-86.2%), and recreation (-81.5%)diminished drastically, much more than any other ecosystem service functions (Table 8).In the study area, therefore, expansion of agricultural land for food production was the main driver causing ESV changes in the area.Ecosystem coefficient of sensitivity (CS) analysis.As shown in Table 9, from 1984 to 2021 the adjusted 50% sensitivity coefficient value of each land use type ranged from 0 to 0.68, all being less than 1, which means the estimated ESVs are inelastic (low sensitive) with respect to the value of modified conservative ecosystem coefficients.Concerning each LULC type, the CS ranged from a low of 0.09-0.2(9.1-20%) for grassland and forestland, to a high of 0.28-0.68(28-68.7%)for water bodies and cultivated land (Table 9).However, the CS value of the water bodies was dramatically reduced from 0.68 (68.2 and 68.7%) in 1984 and 1991 to 0.32 (32%) in 2001 and 0.28 (28) in 2021 (Table 9).Meanwhile, CS for cultivated land exhibits an increment over the observed period .In this study, the barren land and built-up area adjusting value coefficients had no effect on ESV, and hence the sensitivity coefficient was equal to zero.Overall, the findings show that the ESVs' estimation was accurate and practical.

LULC dynamics: Scenarios to the status of dry evergreen Afromontane forest
In the last 37 years, the study areas have undergone significant changes in land use and land cover, largely as a result of the conversion of forests and grassland into agricultural land and built-up area (Figs 2 and 3 and Table 6).According to the LULC analysis of changes in spatial patterns of forest loss over time, forest loss significantly increased in the second phase (1991-2001) compared to the first (1984-1991) and third (2001-2021) phases (Figs 2-6), which is primarily attributable to anthropogenic activities.Similar findings of declining forestland cover have been reported in various parts of the country [e.g., 14, 16-18, 25, 43, 66, 67].These studies concluded that accelerated population growth and subsequent urbanization proceeded to drive uncontrollable resource utilization, and expansion of agricultural lands and settlements at the expense of forestland, which led to the LULC transitions.In addition, some researchers [e.g., 15,17,64,66,68] explained that since there was no government oversight during the interim period that followed the fall of the Derg Government of Ethiopia in 1991, illicit logging activities had occurred throughout the country.Also, this may have contributed to the decrease in forest cover in the study area during second period.However, the percentage of forest loss has gradually decreased, particularly in the third time frame from 2001 to 2021 (Figs 2 and 3).Additionally, NDVI analysis confirmed an increase in the proportion of healthy vegetation during the last two decades (Figs 5 and 6).This may be attributed to Ethiopia's Sustainable Land Management (SLM) and Climate Resilient Green Economy (CRGE) policy initiatives enacted in 2008 and 2011, respectively [34,69].As part of these efforts, the government has been implementing a range of sustainable forest management intervention programs, including promoting soil and water conservation (SWC) measures (e.g., bench terraces, trenches, soil and stone bunds, etc.), planting different trees, shrubs, and herbs, establishing enclosures/exclosures, and discouraging free grazing, hence ensuring a sustainable flow of ecosystem services, as it was also reported in other similar studies [15,65,66].This study also confirmed a lower annual rate of forest loss (1.44%) when compared to other studies conducted in dry evergreen Afromontane forests at various sites as a result of the active community involvement in sustainable forest management activities.For instance, Minta et al. [16], Hishe et al. [67], and Demissie et al. [68]   Our finding indicated that grassland was among the LULC types that changed the most rapidly during the study periods (Fig 3 and Table 6).This finding might be related to the community's rising demand of land for farming, settlements, and grazing, as well as to institutional and policy considerations.For instance, the Amhara Regional State decreed land redeployment in 1991, allocating a large portion of the grazing lands to younger landless peasants for agriculture.There have been reports of large-scale conversion of grassland into agricultural land and settlements throughout the country [16,44,46,68].Gebrehiwot et al. [20] have also reported in the nearby agro-ecology of Ethiopia's Abune Yosef mountain range, a decrease in shrub and grasslands and a commensurate increase in agricultural and built-up areas between 1986 and 2017.Additionally, the progressive expansion of a built-up area in the study landscape may be attributed to the rapid population growth following rapid socioeconomic transformation.According to the Central Statistical Agency [70], the study area had 164358 residents in 1994, but that number alarmingly grew to 197444 in 2007 and is projected to rise by 0.5% annually until 2021.Urbanization will undoubtedly require land for infrastructure and habitation [65], and hence the expansion may be employed at the expense of other land use types, such as natural ecosystems (grasslands, forests, and water bodies).The area's extensive barren land is linked to the parent materials, which are of volcanic and terrestrial sedimentary rock origin [71], as well as to previous incidences of massive erosion and deforestation.This finding is consistent with Gebrehiwot et al. [20], who have reported a comparable large proportion of barren land in the surrounding district of the Abune Yosef mountain range.Similar sharp declines in the water bodies have been seen in a variety of Ethiopian landscapes.For instance, the trend of water bodies showed declines of 19.2% in western Ethiopia [25], 97.5% in the Didessa sub-basin of Ethiopia [64], and 33% in the upper Blue Nile basin of northern Ethiopia [72].The study further showed that a substantial portion of the study landscapes experienced changes in LULC.Hence, the LULC transition matrix in Table 6 reveals the forestland, grassland, and water bodies have been highly converted into agricultural land and settlement areas.This is in line with other studies conducted in Ethiopia [e.g., 16,17,65], which described a marked increase in the proportion of cultivated and built-up areas that could be attributed to satisfy the demand of the expanding population for food production, shelter, and urbanization.Designing and implementing appropriate land use policies in line with the needs of settlement and agricultural land is therefore essential to bringing about unwavering regional sustainability.

Impact of LULC dynamics on ecosystem services and functions
Globally, LULC transitions are changing rapidly, mainly owing to human-induced adversity, such as population growth, urbanization, and agricultural expansion which in turn, shrink the global ecosystem services [73,74], and costs the world over USD 6.3 trillion, equivalent to 8.3% of global GDP in 2016 [75].According to Costanza et al. [8], the global total terrestrial ESV decreased by 28.82%.Our estimation shows that the region's total ESVs decreased by 32.7%, from an estimated value of 80.  7).This estimate is less than both the adopted local value coefficient by nearly two times and the global estimate of Costanza et al. [8].The drastic decrease in total ESVs can be attributed to an expansion of settlements and cultivated land at the expense of natural ecosystems with high value per hectare (e.g., forestland, water bodies, and grassland).Similar losses in total ESV have been reported in various areas of Ethiopia [e.g., 29, 31-33, 65, 72] and even in some other countries like Kenya [58], China [73], Central Asia [74], Nigeria [76], and India [77], who have reported the loss of water bodies and forestland as the most LULC types contributing to the reduction in total ESV.Despite ESVs that showed a significant decline over the entire stipulated period, the proportion of ESV loss attributable to forest loss has gradually decreased over the past two decades (Fig 7 and Table 7).This is probably a result of the government's execution of several sustainable forest management programs, which aim to restore ecosystem health in degraded areas.This is consistent with previous study reports in Ethiopia by Tesfay et al. [65], who have found that the establishment of area exclosures on degraded land has a positive effect on the ESV, with an increasing trend from 2.37 million USD in 2008 to 4.23 million USD in 2017.Importantly, the large increase in cultivated land may result in the loss of natural ecosystem services, despite appearing to be economically profitable [74].Although several interventions and management options have been implemented in the area, the influence of cultivated land has been considerable, leading to decreases in the ESVs of forestland and grassland of 12.2% and 34.5%, respectively.Similar losses of annual ESVs due to agricultural land expansion at the expense of forestland and grasslands have been reported by Kindu et al. [29] in the Munessa-Shashemene landscape of Ethiopian highlands, Tolessa et al. [30] in the Toke Kutaye district, central highlands of Ethiopia, Solomon et al. [19] in the dry Afromontane forest of Northern Ethiopia, and Anley et al. [32] in the Rib watershed, Upper Blue Nile Basin, Ethiopia.
In response to individual ecosystem services, the results confirm that most of the ecosystem service functions have incessantly dropped over the study periods.This decline is mostly attributable to LULC changes that converted grassland, forestland, and water bodies into agricultural areas.As a result, the water supply, hydrological regulation, recreational, waste treatment, and disturbance regulation are largely affected during the stipulated periods.This could, in turn, have an impact on the resilience and resistance potential of the area to ongoing human-induced factors, as well as the livelihood of individuals who depend on the land.We compared our results with those of previous studies on the influence of LULC change on ESV in Ethiopia.And, these studies found that, despite continuous increases in the production of food, timber, and housing, human-dominated land-uses (such as urbanization and agricultural demands) have had a detrimental impact on the provision of ES [29-31, 33, 72].This suggests understanding the drivers of changes and calling for intervention to maintain natural habitats (forestland, water bodies, and grassland) for climate regulation, hydrological regulation, water supply, water treatment, and recreation roles.In the observed periods, the inelastic ESV indicates that the results are reliable, as recommended by Li et al. [62] and Kindu et al. [29].The highest CS value is found in water bodies because the associated value coefficient is relatively higher, while the highest value is found in cultivated land since there is more cropland in the study region than there is in other land use types.

Limitations of the study
In this study, historical LULC changes and ESVs estimation were examined using more than three decades of Landsat imageries.However, remote sensing imageries do have some limitations.For instance, the data between 2003 and 2013 were stripped; it was extremely challenging to assess the LULC trajectories even after de-striping.We used the benefit transfer value approach to estimate the ESVs, and employed a modified conservative value coefficient adopted from Kindu et al. [29], which works based on the assumption that the value of ecosystem services would be uniform across all biomes over time.However, most ecosystems are diverse in terms of their services and functions [8,29,30].Additionally, Kindu et al. [29] noted that the modified conservative value coefficient is uncertain and the ESV in settlements and barren land was not considered in this study because there was no suitable equivalent factor.However, the large-scale expansion of the built-up area in this study had an impact on the natural environments, which has the greatest ESV value coefficient.Despite these drawbacks, the assessment of spatiotemporal LULC trajectories and the estimation of ESVs in response to LULC trajectories are important in providing information for regional level decision-making and in improving the accuracy of these estimates.

Conclusions
Our findings provide evidence that long-term human influence on the environment has resulted in a major change in LULCs in north-eastern highlands of Ethiopia during the past thirty seven years.Because of this, the built-up and agricultural areas have undergone the most changes at the expense of natural ecosystems (e.g., grassland, forestland, and water bodies).The large-scale LULC conversion may have been spurred by the accelerating socioeconomic development and rising demand for land for housing and agricultural/food production.Surprisingly, the rate of forest loss gradually decreased throughout the third study period (2001-2021), which may be the result of government efforts to manage forests sustainably.According to our research, agricultural ecosystem services increased by about 2.2 US$ million between 1984 and 2021, with a 29.7% increment per annum.Together, the total areas of forestland and water bodies dropped by -5.7 US$ million (53.8%) and -47.7 US$ million (86.6%), respectively.This substantial increase in agricultural land and decline in the natural environment may result in a decline in ecosystem services and functions (e.g., water regulation and climate regulation).Thus, policy makers must take into account the major drivers associated with the broad-scale expansion of cultivated land and built-up area and the subsequent diminution of natural ecosystems to deprive and curb further ecosystem degradation and concomitant losses of ecosystem services.As far as we are aware, evaluation of remote sensing-based analyses of land suitability and forest restoration status at a regional level and national-level are quite rare.Our findings thus pave the way for additional research into the prediction of future LULC dynamics nexus ESV using remote sensing observations, modelling land suitability analysis and forest restoration status, and the development of more accurate and trustworthy value coefficients.Moreover, using this baseline data, further studies are suggested that focus on the evaluation of the economic repercussions brought on by the LULC transition at the regional level in particular and in other geographical settings throughout the country.

Fig 3 .
Fig 3.The proportion of LULC types and the area changed in hectares between 1984 and 2021.The positive changes suggest an increase/gain and the negative values indicate a decrease/loss in the status of land use types.https://doi.org/10.1371/journal.pone.0289962.g003

Fig 4 .
Fig 4. LULC changes indicating area change in percent shares between 1984 and 2021.Note: CL, Cultivated Land; BL, Barren Land; BA, Built-up Area; FL, Forest Land; GL, Grass Land; WB, Water Bodies.The positive changes suggest an increase/gain and the negative values a decrease/loss in the status of land use types.https://doi.org/10.1371/journal.pone.0289962.g004 The values in each of the cells signify the amount of land converted from one land cover type to another; the bold diagonal elements indicate the area of each LULC class that remained unchanged between 1984 and 2021.The Net change (Nc) is equal to the gain minus loss and the net change-to-persistence ratio (Np) is the net change divided by diagonals of each class that did not change[30,64,65].https://doi.org/10.1371/journal.pone.0289962.t00634.96% (34785ha) and 34.8% (34507.4) of the area in 1984 and 1991, respectively.Nevertheless, the area had drastically decreased to 27.6% (27380.4ha)by 2001, and then the vegetation cover had increased to 32.0% (31780.6ha)after two decades (Figs5 and 6).Meanwhile, very low vegetation (NDVI: -0.03 to 0.15) indicates barren land and settlement, its percentage grew from 43.3% (43030.4ha) in 1984 to 46.9% (46517.1ha),56.9% (56517.1ha),and 49.2% (48809.7ha) in 1991, 2001, and 2021, respectively.Consistently, the non-vegetation area (NDVI: -1 to -0.03) increased from 1.5% (1564.1ha) in 1984 to 6.1% (1629.9ha)and 4.4% (4369ha) in 2001 and 2021, respectively (Figs5 and 6).

Table 4 . Ecosystem service values of functions of each LULC type based on the Economic Ecosystem and Biodiversity valuation database [29]. Functional types of ecosystem services Accounting index LULC types ecosystem service values (US$/ha/yr) Forest DAFs Cultivated land Built-up Area Grassland Barren lands Water bodies
Notes: WS, Water Supply; FP, Food Production; RM, Raw Materials; GR, Genetic Resources; Gr, Gas regulation; Cr, Climate regulation; Dr, Disturbance regulation; Hr, Hydrological regulation; EC, Erosion Control; BC, Biological Control; WT, Waste Treatment; NS, Nutrient Cycling; SF, Soil Formation; Ha, Habitats; Po, Pollination; Re, Recreation; Cu, Cultural https://doi.org/10.1371/journal.pone.0289962.t004